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Abstract Magneto-atmospheres with Alfven speed [a] that increases monotonically 
with height are often used to model the solar atmosphere, at least out to several 
solar radii. A common example involves uniform vertical or inclined magnetic field 
in an isothermal atmosphere, for which the Alfven speed is exponential. We address 
the issue of internal reflection in such atmospheres, both for time-harmonic and 
for transient waves. It is found that a mathematical boundary condition may be 
devised that corresponds to perfect absorption at infinity, and, using this, that many 
atmospheres where a(x) is analytic and unbounded present no internal reflection 
of harmonic Alfven waves. However, except for certain special cases, such solutions 
are accompanied by a wake, which may be thought of as a kind of reflection. For 
the initial-value problem where a harmonic source is suddenly switched on (and 
optionally off), there is also an associated transient that normally decays with time 
as or C(t _1 In f), depending on the phase of the driver. Unlike the steady-state 

harmonic solutions, the transient does reflect weakly. Alfven waves in the solar corona 
driven by a finite-duration train of p-modes are expected to leave such transients. 

Keywords: Waves, Magnetohydrodynamic; Waves, Alfven 



1. Introduction 

Alfven waves in a stratified atmosphere are governed by the standard linear wave 
equation 

where £ is the plasma displacement (which is transverse to both the magnetic field 
and the direction of inhomogeneity x), and a = a(x) = \B x \j ^fp is the Alfven speed, 
or more properly the Alfven velocity component in the i-direction (the magnetic 
permeability is scaled to unity throughout). In large-scale open field regions of the 
Sun's atmosphere, the Alfven speed increases monotonically with height due to the 
decreasing density [p] , with geometric diminution of magnetic- field strength playing a 
secondary role. Eventually, a reaches a maximum at several Rq and decreases with 
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distance thereafter (An et al., 1990 



become important and Equation (II 



(Velli, 19931 



. By this stage though, solar-wind flows have 
must be modified to take these into account 



However, our focus in this article lies below such heights, and with wave periods of 
minutes rather than the hours common in the heliosphere. We simply ask, What is the 
nature of Alfvenic solutions of the simple Equation |7p? We are all very familiar with 
the basic wave equation, treated at length in every text book on partial differential 
equations (PDEs). However, the nonuniform Alfven speed introduces some features 
that are perhaps less well known, and that have significance for the nature of Alfvenic 



oscillations now seen in the solar corona (Tomczyk et al., 2007 Mcintosh et al, 2011 ). 



We shall address a range of analytic Alfven speed profiles, but the most basic is 



ai(x) = ao e 



x/2h 



(2) 



which pertains to a uniform magnetic field and an isothermal density stratification 
(scale height h). This model has been much used as a fundamental representation 
of Alfven waves in stellar atmospheres at least since Ferraro (1954). Ferraro and 



Plumpton (19581 noted the exact solution 



A 1 Jo ( 2 ^e-*/ 2 " 
a 



+ A 2 Y Q 2 



uh 
— f 
a 



-x/2h 



(3) 



for a wave of single frequency [ui] in terms of Bessel functions of the first and second 
kind of order zero, with Ai and A-i arbitrary constants. They then dropped the Yq 
solution on the grounds that the velocity perturbation does not vanish as x — V oo (i.e. 
as p — > 0). We term this the regularity boundary condition. Many subsequent studies 
have adopted the same approach ( |An et al., 1989 for example). This is problematic 
though; it imposes a perfectly reflecting boundary at infinity thereby setting up 
a standing wave J Q (2e~ x / 2h uih/ao) , which may not be desirable or realistic. The 
regularity boundary condition is unnecessary from an energy point of view too. 
Despite the Y solution being unbounded, 0{x) in fact, the kinetic-energy density 
[|pa; 2 |£| 2 ] vanishes at infinity, and the magnetic energy density \\b\ 2 is finite there, 
where b is the magnetic field perturbationj^The £ = 0(x) behaviour is just the linear 
flapping of rigid field lines to be expected physically as a — > oo (see the animation 



attached to Cally and Hansen 2011 ) 



Although mathematically we are at liberty to impose the regularity boundary con- 
dition, it is often more convenient in theoretical modelling to allow waves to escape 
at the top so as not to confuse matters with downward travelling waves reflected from 
an unphysical "infinity". The boundary condition at infinity is particularly important 
in the exponential model since the Alfven travel time [r = 2h/a(x)] from any point x 
to infinity is finite. To place this in a solar coronal context, assuming a 2 G magnetic 
field, a base density of 10 -12 kgm -3 , and a density scale height of 20 Mm, the travel 



A quadratic wave-energy equation is easily constructed from the linearized momentum and 
induction equations: dE/dt + dT/dx = 0, where £ = \pv 2 + \b 2 is the energy density, 
T = — B bv is the wave-energy flux, v = d!;/dt is the plasma velocity, and b = — B d!;/dx is the 
magnetic field perturbation. Energy fluxes may be attributed to solutions of the wave equation 
using the formula for J 7 , and hence reflection coefficients may be calculated if these solutions 
can be split into upward- and downward-propagating parts. If £ and b are being modelled as 
complex quantities, then we should write £ = |p|«| 2 + ||&| 2 and T = —B Re(f>i>*). 
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time to infinity in an exponential atmosphere is around 450 seconds, comparable to 
the period of the waves that we are considering. 

Despite Alfven waves undoubtedly reflecting from sharp features such as the 



chromosphere-corona transition region (Cranmer and van Ballegooijen, 2005 Hansen 



and Cally, 2012 ), or returning along closed loops, there is reason to believe that waves 



in the open-field corona are preferentially outgoing. This is partly because some 
fraction of their energy irrevocably escapes into the solar wind, and partly because 
of damping by such mechanisms as turbulent cascade to the ion-cyclotron scale (see 



for example Marsch, Vocks, and Tu 



2003 



Hollweg 2006 1, neither of which we choose 



to model here. So, a radiation condition is a way of avoiding such complications when 
exploring propagation far below where they are situated. 

The wave-energy flux associated with the general solution ^ is 



UjB 2 

ihw 
ujB 2 

h 7T 



\A 1 + iA 2 \ 2 



(4) 



where we define Si = |(^4i — iA 2 ) and S 2 = |(^4i + iA 2 ). This shows that the 
alternate but equivalent representation 



uh 
2 — e" 
a 



x/2h 



ujh 
2 — e" 
ao 



x/2h 



(5) 



separates the solution into downgoing and upgoing parts respectively. The solution 
with Hj = 0, 



£f - if. 



(2) 



ujh 
2 — e" 

a 



x/2h 



(6) 



was used by |Schwartz, Cally, and Bel ( 1984 1 to represent an outward-travelling wave 
(the harmonic temporal dependence is omitted from now on, but is implied). This 
Hankel function (Bessel function of the third kind) is well known to asymptotically 

(2) / 

reduce to a complex exponential Hq (r) ~ •v/2/(7rr) exp[— i(r — tt/4)] for large 
argument r ->■ oo (x -> -oo) ( |NIST, 2010] formulae 10.17.5-6) ( |Hollweg| ( fL978 ) 
also uses a Hankel-function solution to represent a radiation boundary condition, 
but assumes that this is valid only "sufficiently far from the Sun that the eikonal 
approximation is valid and no more wave reflections are expected".) 



Hansen and Cally (2012 1 point out that the Alfven wave equation with exponential 



Alfven speed ai(x) is isomorphic to the uniform axisymmetric two-dimensional (2D) 
wave equation with unit wave speed 



9^ 
dt 2 



id_ ( di 

r dr \ dr 



(7) 



under the transformation r = 2h/a. Here the axis r = corresponds to x = +oo in 
the Alfven wave problem. The textbook solution J (ujr) = | [//^'(cjr) + '(ur)] 
effectively imposes perfect reflection on the axis, resulting in a s tanding wave. Placing 
a harmonic source there instead yields the Hq (ojr) solution (Courant and Hilbert 
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1962 Ch. Ill, Section 3.2; Whitham 1974 Section 7.4) whilst the time reverse corre 



sponds to a perfect absorber and yields Hq '(cur). There is nothing reflective in this 
system. But neither are the waves represented by simple d'Alembert-like solutions. 



As is well known and understood (Hadamard, 1923 Courant and Hilbert, 1962 



pp. 208-210), solutions of the wave equation for uniform media in even-dimensional 
spaces do not obey Huygens' principle. Any propagating disturbance trails a wake, or 



"reverberation"; to quote John ( 1982 I, "Disturbances propagate with finite speed but 



after having reached a point never die out completely in a finite time at that point, 
like the waves arising from a stone dropped into water". But even more restrictedly, 
only in one and three spatial dimensions do spherical waves propagate "relatively 



undistorted" (Courant and Hilbert, 1962 Ch. VI, Section 18). The 2D case of interest 



here necessarily exhibits a wake. 



As noted by Hansen and Cally (2012 I 



(1,2) 



(r) 



1 e ± 1U r 2i 

= du T — 

2 ^ TV 



Wo v 7 ! 



-r sinh ', 



dr 



(8) 



Clearly, a + corresponds to propagation away from r = (backward in x), a_ 
represents propagation towards r = 0, and /3 is the reverberation. We see that the 
propagation part of £(r), i.e. a±(ur), consists of a superposition of all wavenumber 
components between and oj. Of course, /J may also be represented as a Fourier 
integral 



/3(r) = 



sin ur du 



+ 



cos ur du 



Vu* 



1 



(9) 



showing that the reverberation can be represented as a superposition of standing 
waves, in agreement with the characteristics-based findings of |Hollweg and Isenb erg 
(2007). In terms of Bessel and Struve functions, a±(r) = Jo(r) ±iH (r) and /3(r) = 
H (r) — Y (r); note that a± is regular at r 
10). For small r, the reverberation 
from Y , /3 • 



(Hansen and Cally, 2012 Figure 



3.1 



and 



3.2 1 in the 



r, tne reverberation [0\ exhibits a logarithmic si ngu larit y in herited 
-27T -1 lnr. This will be seen again later (Sections 
asymptotic form of the transients. 

The fact that Alfven waves in the inner corona appear to be preferentially out- 

(2) 

going suggests that Hq {cor) is a better model of solar atmospheric Alfven wave 
propagation than is Jo(wr), and an indication that strong absorption occurs before 
the waves can reflect. In any case, we are mathematically at liberty to impose a 
radiation condition at large x, and it is certainly convenient to do so in theoretical 
investigations of outward wave propagation. This article, in part, addresses the con- 
sequences of this choice, both for the steady-state harmonic wave problem and for 
the initial-value problem. 

Indeed, it is very common to seek to impose a radiation condition above an 
exponential or similar atmosphere by simply appending a uniform plasma above 



some (large) height (Hartree, 19291. However, it will be shown that this is not a 
good choice in general because the discontinuity in Alfven speed gradient is itself 
highly reflective. The artifice therefore decides the issue rather than illuminates it. 

Another interesting suggestion for how to avoid the difficulties presented by a 
finite travel time to infinity and perfect reflection there is to retain the displacement 
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current in the wave equations (Leroy, 1983 Tsap, Stepanov, and Kopylova, 2009). 
Then the Alfven speed is limited by the speed of light, and the Alfven wave ultimately 



couples to an outgoing electromagnetic wave. (It should be noted though that Leroy 



did not mean to imply that this process actually occurs in the solar atmosphere - see 
the last sentence of his article - the device was merely introduced as a mathematical 



way around the troublesome infinity.) Unfortunately, we shall see in Section 2.3 that 
this model is almost totally reflective for any realistic frequency and density scale 
height, and so does not fulfil our requirement for optimal transmission. 

Reflection of harmonic waves will be addressed from several standpoints and for 
a variety of atmospheres in Section [2] showing that the reflectivity of an atmosphere 
is governed by its smoothness, or lack thereof. Discontinuity in any derivative of 
the Alfven speed results in a reflection coefficient that depends algebraically on the 
wavenumber. On the other hand, an infinitely smooth ( < ff°°) function suffers only 
exponentially small reflection at worst. Then we address the initial- value problem, 
in the exponential atmosphere, with specific focus on the decay of transients, both 
for the radiation and regularity boundary conditions. Finally, we draw some general 
conclusions about the relationship between reflectivity, wakes, and transients and 
how they are determined by the form of a(x). 



2. Reflection 

2.1. Power Law Alfven Profiles and the Prevalence of Wakes 

One might expect that any nonuniform wave-speed profile gives rise to reflection. 
This is not the case though. For example, it is known that the profile a = x 2 admits 
relatively undistorted exact solutions of the form £ = xU(t ± x~ v ) that represent 



unidirectional propagation ( Didenkulova, Pelinovsky, and Soomere, 2008| . 

For other power laws a = aoX < > n ~ 1 - >/< - n ~ 2> (n ^ 2, with x > if required) the wave 
equation transforms to 

dt 2 dr 2 r dr 

under the change of variables r = {n — 2) aQ 1 x~ 1 ^ n ~ 2 \ The profile a(x) is monotonic 
increasing if n > 2 or n < 1, with the former case yielding finite Alfven travel time to 
infinity, J a(x) _1 dx < oo. For n $C 1 the Alfven exponent (n— 1) /(n—2) 6 [0, 1) and 
the travel time to x = +oo is infinite, so the issue of reflection from infinity does not 
arise in any initial-value calculation. For positive integer n (although we have no need 



to restrict it thus) Equation ( 10 1 is the spherical wave equation in n dimensions. In 
one dimension the well-known d'Alembert solution is £(x,t) = U(t±r) for arbitrary 
shape function U (this is strictly undistorted) , and the analogue for three dimensions 
is £{r,t) = r~ l U(t ± r) (relatively undistorted). The n = 3 case corresponds to the 
a = x 2 profile mentioned above. 



Exact harmonic solutions of Equation ( 10 I may be represented in terms of Hankel 
functions, 

t = r»H { ^\ur)e-^\ (11) 

where /i = 1 — (n/2). These solutions reduce to elementary functions without branch 
points only when n is an odd integer (spherical Bessel functions). For n = 1 (a = 1, 
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arbitrarily scaling a to unity), £ = e lu ( ±r 



, from which the D'Alembert 



solutions clearly may be constructed by Fourier composition. For n = 3 



2 ). 



i±x~ 1 -t) 



which again obviously leads to £ = x U(t ± x 



), £ = r 3 (wr =F e J " (t±r ^ with r = 3x 1 / 3 , showing that 
= r~ 3 U±(r ± i) — r~ 2 U±(r ± i) for arbitrary differentiable 



£ = r -l e iw(±r-t) = xe k 

for arbitrary U. 

For ?i = 5 (a = x 4 / 3 
solutions of the form £ 
C/± may be synthesized. Initial displacement specifies U+ + U- and initial velocity 
gives U+ — U-. If £(r, 0) and £ t (r, 0) have joint compact support, £(r, t) will continue 
to do so, with the interval on which £ ^ simply translating at unit speed in each 
direction. No wake is left, although the shape of the pulse changes. Higher odd-n 
cases behave similarly, with progressively higher derivatives of U involved. In fact, 
£ = (d 2 /dt 2 — V 2 )( n_1 '/ 2 (/>(t ± r) is a solution for arbitrary 0, where V 2 is the n- 



dimensional spherical Laplacian, i.e., the right-hand side of Equation ( 1 1 (Courant 



and Hilbert, 1962 Ch. VI, Section 13.4). 



For all other real powers, and for the exponential profile that corresponds to 



?i = 2 in Equation (10 I, radiating solutions contain logarithmic or fractional power 
singularities at r = 0. These branch points give rise to continuous spectra in the 
initial value problem, and hence to transients that decay algebraically with time. 
Simple "relatively undistorted" translating solutions exist only for n = 1 and 3, and 
compact support is simply translated for other odd n. Otherwise, propagating wave 
packets trail wakes. This is the rule rather than the exception in inhomogeneous 
media. But do we get reflection other than that implicit in the wake? 



2.2. WKB Perspective 



Within the WKB approximation ( |Bender and Orszag, 1978 Gough, 20071, one seeks 



asymptotic solutions of a £ X x + w 2 £ = that take the form 



£ = exp 



u;- m S n (x) 



(12) 



Truncating at m = 1 (the physical optics approximation), the general solution for 
arbitrary a(x) is given by 



£ - a(x) 1/2 [Ae llp(x) + Be~ iv(x) } 



uj — > co, i.e. k — > co, 



(13) 



formally valid for k ^> |dlnfc/da;|, where k(x) = ui/a(x), fix) = J x k(x) dx = —iSo, 
and A and B are arbitrary constants representing the amplitudes of the upward and 
downward waves respectively. The two waves are completely decoupled, meaning 
that there is no reflection. Coupling is not recovered by going to higher order in to. 



The reason for this is that the original assumption ( 12 1 does not admit such coupled 
solutions. 

However, this all breaks down where k~ l |dln k/dx\ 1, and particularly at points 
of discontinuity in a or a'. For discontinuous Alfven speed, the wave-energy reflection 
coefficient results from the m = term and is given by 1Z = (a + — a_) 2 /(a + + a_) 2 
to leading order, independent of frequency, where a± are the Alfven speeds on either 
side of the discontinuity. For discontinuous Alfven speed slope the effect enters at 
m = 1 and yields TZ = (Aa') 2 /(16tJ 2 + (Aa') 2 ), where Ao' = a' + - a'_ is th e jump 
in slope. Notice that reflection vanishes as u) — > oo in this case. Velli (19931 makes 
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the interesting point that multiple such discontinuities in scale height can set up 
resonances because they create (leaky) cavities. Discontinuities in higher order [d > 
1] derivatives are invisible at the physical optics level, appearing only for higher 
m = d and yielding 0(u~ 2d ) reflection coefficients. We therefore expect that 1Z 
is exponential rather than algebraic for < <f 00 functions a(x) (i.e. functions that are 
continuously differentiable to all orders). This is confirmed in the examples below. 

2.3. Asymptotically Flat Profiles 

Consider several cases where the exponential Alfven profile for a 2 is capped at a 
maximum level a 2 , = ag/e as x — > oo. 

• Sharply Capped Exponential: For 



<Jq e x l h x < —/line, 



ag/e 



x > —hhi e. 



(14) 



the reflection coefficient is 



Ki = 



H^ 2) (2k u h) + iH^>(2k u h) 



(2)/ 



H^ ) (2k u h) + \H\ L> {2k u h) 



(i). 



1 



MB.h 2 



(15) 



as k u —¥ oo, where k u = uj/a u , in agreement with the matched WKB result 
above. 

Smoothly Capped Exponential: With a 2 = a\ = a 2 , e X/ '' l /(l + ee 1 / 1 ), the 
harmonic wave equation can be expressed in the form of a Bessel equation, 



(16) 



where k = u>/a Q = e~ 1 / 2 k u and r = 2k he~ x / 2h . The solution satisfying the 
radiation condition at x — > +oo is easily identified using J v (z) ~ (z/2) v /T(\+v) 
as z — > 0: 



£ = 2J_ 2 ifc u /i(r) since it is 0(e tkuX ) as x — > +oo 

(r) 

-irkuh e i(r-it/i) _|_ e 7r k u h g-i(r-vr/4) i ag r 



2J-2ik 


A 






m 






e~ 


V 7r r | 





(2) 

—2 i k u h V' 1 (17) 



The reflection coefficient is therefore the squared ratio of the coefficients of the 
complex exponentials: = e -47rfc "' 1 . As expected since a(x) is analytic, reflec- 
tion is exponentially small as k u — > 00. Equation (161 has a regular singular 
point at r — (x — 00) with roots of the indicial equation /1 = ±2ik u h. 
Alternate Smoothly Capped Exponential: With a gentler profile a 2 = a\ = 
a 2 e x/h _|_ £ i/2 e x/2hj2 tne wave equation takes the form 



d ( d£ 
r— r 



dr \ dr 



+ (r + k u hf i = , (18) 
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Figure 1. Transmission coefficients against k u h for the three asymptotically flat profiles. Ti 
(full curve); T3 (dashed); 74 (dotted). All three transmission coefficients are asymptotic to 
4izk u h as k u — > 0. 



with r = 2k he x / 2h again. The radiation solution is given in terms of 
confluent hypergeometric function, 



xFx [ l;l-4ik u h;2ir 



(19) 



Analysis of the large r asymptotics of the 1-F1 function (NIST, 2010 formula 
13.7.2, where M(a, b, z) = i-Fi(a; b; 2) /]?(&) )reveals that the reflection coefficient 
is IZ4 = e~ i7rkuh sech(4nk u h) . Once again, 1Z decreases exponentially with in- 
creasing k u , as expected. In this case, r = is again a regular singular point, 
but now the roots of the indicial equation are p, = ±i k u h. 

The three transmission coefficients T = 1 — 1Z are plotted against dimensionless 
wavenumber [k u h] in Figure [l] illustrating the fact that 'if 00 profiles remain essen- 
tially fully transparent to much smaller wavenumbers. Lest the reader think that the 
differences are entirely due to the maximum values of the WKB validity parameter 
A; -1 |d In A;/dx| , i.e. (2k u h)~ 1 , (3^k u h) , and (8fc„/i) -1 respectively for the three 
cases, this does not fully account for the discrepancies. 

The three capped models are all totally reflective in the limits a u — v 00 and 
uj — > both of which correspond to k u —5- 0, so they represent poor devices for 



allowing radiation to escape at the top. This approach led Bel and Leroy (1981) 



to discount Alfven losses in sunspots. As suggested by the figure, a sequence of < tf°° 
functions, more and more closely approximating the sharply capped exponential, will 



of course return the transmission coefficient (15 1 in the limit. 

Oq/c 2 , where c is the speed of light, the smoothly capped exponential 



With 



case corresponds to that of Equation (9) in Leroy (19831, where the displacement 
current is retained in the Alfven wave equation and the wave travel speed is limited 
by c. The corresponding reflection coefficient [TZ = e _47rw,i / c ] is of course near total 
in any conceivable realistic case, since coh/c <S 1. The large outward fluxes found by 
Leroy are an artefact of his prescribed driver £t(x, t) = V e~ lujt at the base x = and 
the near-impermeability of the top setting up a weakly leaky resonant cavity. Then 
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£ = i Vu 1 J-2iuh/c{rQ e x/2h )/J-2iuh/c{ro) where r = 2ujh/a , with the near-zeros 
of the (complex) denominator defining the resonances. The resonant cavity issue will 
arise again in Section |3.3| 

2.4. General Method to Impose Radiation Boundary Condition at Infinity and Why 
the Exponential Atmosphere is Transparent 



Asymptotically flat profiles such as those considered above are ideally suited to the 
calculation of overall reflection coefficients, no matter how steep and non-WKB they 
may become at intermediate x. However, this does not mean that we can determine 
the reflectivity of a section of a non-WKB a by abruptly sandwiching it between 
two uniform or WKB regions. Even a segment of the clearly non-reflective a = x 2 
bookended continuously by two uniform sections will return non-zero reflection, but 
this is entirely due to discontinuities in a'. 

For unbounded Alfven profiles such as ai(x) though, there is no WKB region as 
x — > +oo in which inward and outward waves can easily be identified. However, the 
change of spatial variable to r = e~ x / 2h or similar that maps x = +oo to r = and 
-oo < i < +oo to r > presents a convenient mathematical device for allowing 
complete absorption of a harmonic wave at x = +oo. Equations such as |7]) and 
(10 1, which are symmetric in r, suggest that waves be allowed to propagate through 
to negative r, thereby removing energy from the physical system. However, since 
the solutions generally contain a branch point at r = and a branch cut along 
the negative r-axis, it is necessary to treat r as a complex variable and address 
the issue of whether to adopt the solution above or below the cut. The solutions 
([6} and (11) must in fact be continued below the cut to obtain the flow-through 
behaviour at r = 0, based on th e analytic con tinuation formula for Hankel functions, 
Hj?\ze- i7r ) = -e i>iw H^\z) |NIST, 201o] formulae 10.11.5)0 This represents a 



mirror image outgoing wave as r — > — oo, and is exactly the required energy sink. Both 
H^\ze~ 17T ) (below cut) and H^ \ze ln ) (above cut) split into linear combinations 



of H^\z) and (z) functions |NIST, 2010| 10.11.8 and 10.11.4), breaking the 



symmetry and precluding a similar construction for solutions such as ( 17 1 that arise 
for complex Frobenius index. 

These considerations indicate that it is possible to mathematically prescribe a 
boundary condition that perfectly absorbs harmonic waves at x = +oo. For Alfven 
profiles with wave equation isomorphic to Bessel's equation with real order, the 
appropriate purely outgoing solution is simply the one involving (uir), i.e., 

(2) / . (2) 

Hq (ujt) for the exponential profile and r^H^^(ujr) for the power-law profiles. In 
neither case is there any reflection. Essentially, this is a consequence of analyticity 
and symmetry in the complex plane. 

For an initial-value problem though, the transient may not be perfectly absorbed. 
This is investigated next for the exponential profile. 



2 Alternatively, keep r real but insert a weak frictional force [— ■yd^/dt] on the right-hand side 
of Equation u7p , with < 7 <C to. Then the upgoing solution is (uiryjl + i 7/01). This 
places the argument of the Hankel function below the cut for negative r. 
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3. Exponential Alfven Profile Initial Value Problem 

We now examine two initial value problems in the exponential atmosphere using 
the method of Laplace transformation. First we consider a harmonic point source 
at r = R that switches on at time t = 0, with either the outgoing radiation or 
the regularity boundary condition applied at r = (x = +oo). Then we allow the 
point-source driver to switch off again at t = T. We are interested in the temporal 
dependence of the transient, and whether it is related to the reverberation (wake) 
seen in the steady state solution Equation ([6}. 

Our solution of the initial-value problem for Alfven waves differs from that of 



An et al. ( 1989 1 in several ways: they solve the wave equation numerically rather 
than via Laplace transformation; they impose a base velocity whereas we impose 
a point force (acceleration); and most importantly they specify a totally reflective 
boundary condition at x = +oo. Their solutions therefore reflect from infinity to set 
up a standing wave in finite time. We address this model too, but only to compare 
with the radiation boundary condition results. 

3.1. Point Source Harmonic Switch On 

Consider the Alfven wave equation with a point driver at r = R, 

w = l !(-af) +>*-'m«r-* m 

where again r = 2h/a, F is an as yet arbitrary function of time that switches on 
at t = 0, and a convenient normalization has been adopted. Laplace transforming 
Equation |7]l with initial conditions £(r, 0) = 0, £t(r, 0) = results in 

p^+ 1 -^+ [jJ 2 u=-8R- 1 f(u)5(r-R), (21) 
or^ r or 

where u(r,u) = /" e 1 w *£(r, t) dt = C{£(r,t)} and /(u>) = C{F(t)}. Imposing out- 
going radiation boundary conditions at both r = and oo, the Laplace transformed 
solution may be constructed as a Green's function, 

u(r, u) = 2 7T i f(u) # (1) K> ) H[? ] {uj T< ) , (22) 

where r< = min(r, R) and r> = max(r, R). The complex inversion formula then 
returns 

/oo+ie 
e-^'/H H^(oj r> ) H^icjr^du , (23) 
-oo+ie 

where e is a positive constant sufficient to place the contour above any singularities 
in the complex tj-plane. The subscript "rad" indicates that the radiation condition 
is applied at r = 00 



NIST 



All formulae used in this section relating to Bessel functions may be found in 
particularly the asymptotic formulae 10.17.5-6 and the analytic continuations 10.11 
10.11.7-8. 



(20101, 



5 and 
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Figure 2. Top row: Stack plots of Re(£) against position r and time t = 0.027T, 0.04/7T, . . . , 0.6-7T 
for the case R = 1 and Q = 25. Left: £ rat i; right: £rcg- Notice how in the left panel the 
phase propagates uniformly to the left throughout, corresponding to the outgoing radiation 
condition, whereas in the right frame a standing wave is progressively set up after the wave 
reaches r = 0. The diagonal black lines in the right panel are t = R ± r; reflection is evident 
only for t > R + r. Bottom left: Transient £ t r for the case of the top-left panel. Bottom right: 
Re(ftr(r, t)) for r = 0.002 (full curve), 0.1 (dashed), 0.2 (dotted), and 0.4 (chained) against time 
[t], showing the slow mo not onic decay of the transient. The long-dashed red curve represents 
the asymptotic formula \27t to leading order i -1 Int for the r = 0.002 case. 



Specializing now to the case of a monochromatic harmonic driver F(t) = e~' nt 
with Q, > 0, we have f(u) = i/(w — O), and only require e > 0. With this choice 
of F(t), Re(£) corresponds to driver — 8i? -1 cosf2i and Im(£) to — 8i? — 1 sin fit. The 
integral in Equation ( 23 I may be evaluated by completing the Bromwich contour in 
the upper half-plane (UHP) if r— R\ > t, and in the lower half-plane if r— R\ < t. The 
former yields ^ ra d = since there are no singularities in the UHR This is expected 
as the signal from the source travels at unit speed in r-space. The \r — R\ < t 
case picks up contributions from the simple pole at u = Q, and from the branch 
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cut conventionally lying along the negative real axis for Hankel functions. The pole 
produces the steady state and the branch cut contributes the transient: 



Cmd(r, t) 



oo+ie -iuit 

ffW(wr>)^ J) (wr<)dw 



uj — n 



(24) 



£ te (r,t) + 2 7rie- im ^ 1) (r!r > ) H^(n r< ) \U(t-\R-r\), 



(2), 



where hi is the unit step function and the transient is 



fcr(r,t) = 2 



fl^ax) H^'(xR) - H^'(xr) H^'(xR) 
~J (xr)Y (xR) + J (xR)Y (xr)] dx. 



4i 



d.r 



(25) 



The oscillatory Fourier integral for £ tr given in Equation ( 25 1 is readily evaluated 



numerically using Mathematica's built-in N I integrate function. An alternative non- 
oscillatory formulation, valid only for t > R + r, is 



&r(M) 



8e-y* 
ifl-y 



I (yr)I (yR) + - (K a (yr)I a (yR) + K a (yR)I (yr)) 



dy. 



(26) 



This is numerically simpler and quicker to evaluate, and more convenient for asymp- 



totics. Specifically, we may use Laplace's method (Bender and Orszag, 19781 to show 
that 

8 / At 2 \ 8 / At 2 \ 

+ — (^-2i + , + iln-j iort»R + r, (27) 

correct to 0(t~ 2 ). It is interesting that the asymptotic decay of Re(£ tr ) is monotonic, 
and depends only logarithmically on r (linearly on x), whilst the decay of Im(£tr) 
is independent of r to leading order. The logarithmic dependence on r reflects the 
small- r structure of the reverberation [/3]. 

For comparison, the same problem but with the usual regularity boundary con- 
dition applied at r = rather than the radiation condition yields 



J (wr<) F^ 1) (wr>)dw 







(28) 



A similar calculation was performed by Bogdan and Cally (1997) for Alfven waves 



in polytropic atmospheres with point impulse initial conditions and a regularity 
boundary condition where the density vanishes. Clearly, since 2Jq{z) = H^\z) + 



A(r, t) = C r ad(r, t) - £ reg (r, t) 



ff^(wr<) J ff«(a;r > )d W , (29) 



n 
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Figure 3. Same as the top-left and bottom-right panels of Figure[2] but for Im(£ rat j). In this 
case, as suggested by Equation j27|l, the t~ 1 decay is independent of r to leading order. 



which vanishes for t < R + r since then the contour must be completed in the UHP. 
This confirms that the two solutions at r are identical up until the time that the 
signal in the regular case reflects from the origin and returns to r. This is as it must 



be. Evaluating Equation (29 I by contour integration 



( e 

AM) = -U(t - r - R)\2 J (5J (xr)J (xR) - Y (xr)Y (xR) 

+ i [J (xr)Y {xR) + J (xR)Y {xr)}) dx + 2 Trie"™ (fir) ^(OR)} . (30) 

Figure [2] (top row) illustrates the difference between the two solutions. They are 
indeed identical until t = R, but thereafter the reflected wave in the "regular" case 



progressively sets up a standing wave (An et at, 1989). We also see that, even 



in the "radiation" case, a very slowly decaying transient is set up that spreads to 
progressively larger r (bottom row). 

Figure [3] displays the behaviour of Im(£ m d) and the asymptotic decay of its 



transient, which is now independent of r, as expected from Equation (27 1. 

We mention briefly that the present analysis may be replicated for the smoothly 
capped exponential profile (not presented here), although with added mathematical 
difficulty due to the frequency [uj] appearing in both the argument and the order 
of the Bessel solutions. As expected, this case displays a marked reflection of the 
transient off the transition from exponential to flat Alfven speed, but it also leaves 
a wake that decays as C(£ -1 ). 

3.2. Point Source Harmonic Switch-On/Switch-Off 

Now consider F(t) = e~ int U(T - t), where T > is a finite driving interval. Then 
f(ur) = i(l — e(" _n ' T )/(o; — fl) has only a removable singularity, so there is no residue 
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Figure 4. Left: Stack plot of Re(£y) against position r and time t = 0.027T, 0.047T, . . . , 0.08-7T 
for the unbalanced case R = 1, Q = 25, and T = tt/5. The wake in the "V" at the right 
(R — r + T < t < R + r , T/2 < r < 1) is extremely weak and not discernible on the figure. It 
is non-zero though. Right: Re(^T(r, t)) for r = 0.002 (full curve), 0.1 (dashed), 0.2 (dotted), 
and 0.4 (chained) against time t. Th e red dashed curve depicts the leading order asymptotic 
transient formula from Equation (132b at r = 0.002. 



and therefore no non-trivial steady state solution (of course). Adopting the radiation 
boundary condition again, the analysis is hardly changed and we find 



f (i 



= i( w -nm 



H, 



6ad fat) 



n 



Ida; 



(31) 



This is a simple consequence of the second shifting property of Laplace transforms 



(Spiegel, 19651. The decay of the associated transient is therefore given asymptoti- 



cally by 



vT.tr 



~ mT -iQT) 



-i 7r + In 



4£2 
rR 



it* 



-2i + 7r + iln^— ) +0(t- 3 \nt) for t > R + r. (32) 



The decay is then 0(t~ 2 ) or C(t _2 lni) if £IT = 2n7t for integer n since the two 
individual transients cancel to leading order. We shall refer to this case as "balanced". 
In the unbalanced case, the generic slow 0(t _1 ) or C(t _1 lnt) decay ensues. 

Figure [2] for an unbalanced case shows how the travelling wave signal moves 
leftward and out of the domain, but leaves a rightward-spreading and slowly decaying 
transient long after the driver has been switched off. Figure [5] depicts a balanced case 
where very little wake is left at large t. Nevertheless, although the oscillatory (steady 
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Figure 5. Left: Stack plot of Re(£y) against position r and time t = 0.027T, 0.047T, . . . , 0.08-7T 
for the balanced case R = 1, CI = 50, and T = tt/25. Right: Re(£ T (r,t)) for r = 0.002 
(full curve), 0.1 (dashed), 0.2 (dotted), and 0.4 (chained) a gain st time t. The red dashed 
curve depicts the leading-order asymptotic transient formula | |32[ l at r = 0.002, which is now 
0(t~ 2 lnt). Note also that the amplitude of the transient is far smaller than for the unbalanced 
case of Figure [4] 



state) signal has escaped the domain at r = 0, a non-oscillatory transient is reflected 
in a narrow band in 9 — t space that represents the absence of balance in the initial 
switch-on. It is soon cancelled by the subsequent symmetric switch-off. 

3.3. Prescribed Base Harmonic Displacement 

For completeness, we briefly mention the case where the base (r = R) velocity or 
displacement is prescribed rather than introducing a driving force term as above. 



This is the more common procedure (Hollweg, 19781 Leroy, 1983 Schwartz, Cally, 



and Bel, 1984 An et al, 19891, and leads to probably unrealistic Alfven-frequency 



resonances because the "closed box" has natural frequencies. 

Specifically, letting £,(R,t) = e _int and adopting the regularity boundary condi- 
tion at r = 0, the solution assuming zero initial displacement and velocity is obtained 
by Laplace inversion: 



5M) = 



2tt 



oc+ie 



u> — Q Jq(u>R) 



Jo(jnr/R) 



e -ij„t/R e ij n t/R 
jn ~ + jn + QR 




where the j n (n = 1, 2, . . . ) are the positive zeros of Jq, and < r R. Asymp- 
totically, the j n are nearly uniformly distributed: j n ~ a n + l/(8a„) as n — > oo, 
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Figure 6. Left: Temp oral evolution of Re£ at r = 0.5 for the rigid harmonic driver case 
f2 = R = 1 of Equation \33\ . The first 50 Bessel zeros are retained. The displacement switches 
on at t = ii — r = 0.5 and oscillates wildly due to the many progressively out of phase trapped 
Bessel modes. Right: Re£ vs. r for the same case, at times t = 100 (full) and t = 0.5 (dashed), 
displaying smooth behaviour in space. For this plot, 100 Bessel modes were used. 



where a n = (n — \)tt (NIST, 2010 formula 10.21.19), so simple poles are distributed 
along the whole real line, specifically at u = ±j n /R. In deriving this solution it has 
been assumed that Q ^ j n for any n. The solution for the transient is no longer as 
neatly expressed as before, since there are now no branch points and hence there 
is no continuous spectrum and no integral associated with a branch cut. This is 
because the transient cannot escape the box, either through r = or r = R. Each 
initially excited mode n persists indefinitely as a discrete harmonic of the closed 
cavity < r < R, resulting in a chaotic time evolution but smooth spatial structure 
(Figure [6) . This is a consequence of treating the bottom boundary as rigid for all 
frequencies other than the driving frequency f2. In reality, Alfven waves should be 
able to penetrate into the solar interior due to their large wavenumber in the low 
atmosphere. 

The case with radiation boundary condition at r = is more complex again, as 

(2) 

branch points are involved as well as the complex zeros of Hq . These complications 
are incidental to our main concerns though, and will not be pursued further here. 



4. Summary and Discussion 

We recapitulate our main findings. 

i) Reflection of time-harmonic Alfven waves vanishes in the high frequency (WKB) 
limit provided the Alfven speed profile a(x) is continuous. If a is a t if d ~ 1 function, 
i.e. (d— l)-times continuously differentiable but not more, then the reflection coef- 
ficient 1Z = 0{u~ 2d ) as uj — > oo. If a is a ^f 00 function, then 1Z — > exponentially 
as u) — > oo. 

ii) There exist atmospheres, a = a a;^ n_1 '^ n_2 ' for integer n, where Alfven wave 
propagation is isomorphic under the transformation r = (n — 2) 1 .T -1 /( n_2 ' to 
the n-spherically symmetric solutions of the n-dimensional wave equation with 
unit wave speed in r-space. Only for odd-integer n do disturbances travel without 
trailing a wake. The formulation of the Alfven wave equation as an n-dimensional 
spherical wave equation in a homogeneous medium provides a simple explanation 
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of the Alfven wakes found by Hollweg and Isenberg ( 2007 1 using a characteristics 
and impulse-function approach, and the "continual coupling" inferred by |Leroy| 
( 1980 1 using matrix diagonalization. 

iii) For n not an odd integer, time-harmonic radiation solutions for these power-law 
atmospheres, and the exponential atmosphere which corresponds to n = 2, contain 
branch points and therefore exhibit the characteristics of a continuous spectrum, 
including algebraic time decay of transients. 

iv) Another feature of these solutions is wakes. Wakes may be thought of as a form 
of (continuous) reflection, but are distinct from simple reflection which sends a 
signal back from whence it came. The existence, or otherwise, of branch points 
is a feature of the differential equation, and can be determined by local analysis 
about its singular points without needing to solve the DE. 

v) The atmosphere a = a\ = a e~ x / 2h corresponds to the case n = 2 under the 
transformation r = 2h/a, and so does not admit relatively undistorted solutions. 
The initial-value problem therefore produces a wake that persists even after the 
driver is switched off. However, the steady-state component consists of a unidi- 
rectional travelling part and a stationary reverberation. There is no reflection 
in the standard sense. With the radiation boundary condition applied at r = 
(x = +oo) the steady-state signal freely leaves the domain, but a weak non- 
oscillatory transient remains that decays as £?(£ -1 ) or 0(t -1 lni) (unbalanced 
case), or 0{t~ 2 ) or 0(t~ 2 \nt) (balanced case), depending on the form of the 
driver. This is a sort of slow monotonic relaxation to the equilibrium state. 
Although the steady-state solution suffers no reflection at r = 0, the transient 
does reflect, spreading backwards (ever weakening) at unit speed. 

vi) Partial reflection can be associated with complex Frobenius indices at r = 0. For 
real indices it is possible to construct a radiation boundary condition at x = +oo 
based on the analyticity and symmetry of the wave equation in r that yields no 
reflection of harmonic waves at any frequency. 



5. Conclusions 

The analysis presented here may seem esoteric. Nevertheless, it has practical impli- 
cations for the way that we model wave propagation in the solar atmosphere. There 



are any number of articles (Leer, Holzer, and Fla, 1982 for example) that seek to 



calculate the "intrinsic" reflectivity of coronal Alfven waves by placing a uniform slab 
above an exponential or similar model. However, this says more about the matching 
point than about the underlying atmosphere. The exponential and power-law profiles 
are cases in point. If the atmospheres are allowed to extend unimpeded to infinity 
without truncation, they are entirely transparent. Or, more physically, if an efficient 
and non-reflective wave energy sink is placed high in the atmosphere, there is no 
reflection from the body of the atmosphere either. A simple uniform or WKB slab 
does not represent such a non-reflective sink, although this is not to say that it may 



not be a reasonable representation of the outer corona (Leroy, 1981 1. In the absence of 
any Alfven wave dissipation, such a WKB top would indeed induce strong reflection 
in the underlying atmosphere. But the point is that the underlying atmosphere is not 
necessarily intrinsically reflective (viz., the exponential or power-law Alfven-speed 
profiles). It is the transition to the uniform or WKB top that reflects. 
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This is relevant when performing numerical experiments in truncated model at- 
mospheres. If we choose to, we can postulate a radiation boundary condition at the 
top of our region of interest. In numerical simulations, this is typically done using 
characteristic boundary conditions (Engquist and Majda, 19771 or absorbing layers 
(Berenger, 1994 1. For steady monochromatic waves, the task may be easier, through 
simple matching to a known analytic solution that represents an outgoing wave. For 
the exponential atmosphere, this involves adopting the Hq (ur) Hankel function 
solution (Cally and Goossens, 2008) . The Jo(wr) solution on the other hand is 
appropriate if we want a reflective boundary at infinity. With this in mind, matching 
numerical solutions in a finite domain to a Hankel function or similar radiation 
solution is mathematically well founded and physically interesting. This does not 
imply that the exponential atmosphere really does extend to infinity; it is simply 
an effective mathematical device for imposing a radiation condition at the top of 
a computational domain. A more realistic treatment of solar coronal Alfven waves 
should in fact address the maximum and gradual decline in the Alfven speed beyond 
a few solar radii and the loss of hydrostatic equilibrium that results in the solar wind 
(Velli, 19931, but that is beyond the expository scope of the present analysis. It is 
notable though that | Velli | finds near-total transmission at high frequencies in this 
model, which further supports the use of a radiation boundary condition in simpler 
atmospheres. 

Analysis of the initial-value problem for the exponential atmosphere verifies that 
the steady state Hq (uir) solution does indeed fully depart the physical model 
through x = +oo, though a weak transient remains that decays algebraically in time. 
Since we might expect episodic generation of Alfven waves in the lower atmosphere, 
for example by fast-wave conversion (Hansen and Cally, 20121, the ubiquitous pres- 
ence of such slowly relaxing transients can hardly be avoided. These motions would 
not be identified as waves observationally, as they are not oscillatory in time. 
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